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Abstract 

The Mongolian Steppe is one of the largest remaining grassland ecosystems. Recent studies have reported wide- 
spread decline of vegetation across the steppe and about 70% of this ecosystem is now considered degraded. Among 
the scientific community there has been an active debate about whether the observed degradation is related to cli- 
mate, or over-grazing, or both. Here, we employ a new atmospheric correction and cloud screening algorithm (MAI- 
AC) to investigate trends in satellite observed vegetation phenology. We relate these trends to changes in climate and 
domestic animal populations. A series of harmonic functions is fitted to Moderate Resolution Imaging Spectro- 
radiometer (MODIS) observed phenological curves to quantify seasonal and inter-annual changes in vegetation. Our 
results show a widespread decline (of about 12% on average) in MODIS observed normalized difference vegetation 
index (NDVI) across the country but particularly in the transition zone between grassland and the Gobi desert, where 
recent decline was as much as 40% below the 2002 mean NDVI. While we found considerable regional differences in 
the causes of landscape degradation, about 80% of the decline in NDVI could be attributed to increase in livestock. 
Changes in precipitation were able to explain about 30% of degradation across the country as a whole but up to 50% 
in areas with denser vegetation cover (P < 0.05). Temperature changes, while significant, played only a minor role 
(r 2 = 0.10, P < 0.05). Our results suggest that the cumulative effect of overgrazing is a primary contributor to the 
degradation of the Mongolian steppe and is at least partially responsible for desertification reported in previous studies. 
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The Mongolian Steppe is one of the largest remaining 
grassland ecosystems in the world (Kawamura et al, 
2007) spanning an area of about 1.17 million square 
kilometers (Fernandez-Gimenez & Allen-Diaz, 1999) 
and encompassing roughly 2.6% of the global grassland 
vegetation (Li et al., 2005). While Mongolia's vegetation 
has been described in general (Kalinina, 1974; Lavrenko 
& Karamysheva, 1993) and some plant communities 
have been studied (Pacyna, 1986; Wallis de Vries et al, 
1996), relatively little is known about the ecosystems 
vulnerability to human activity and climate change 
(Sugita et al, 2007). Recent decades have seen significant 
decline in grasslands across the country (Ykhanbai et al., 
2004) and over 70% of the steppe is now considered 
degraded (UNEP, 2002). At the same time, there has 
been a dramatic increase in livestock numbers in Mongo- 
lia (Gong Li et al, 2000); domestic animal populations 
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have almost doubled from 26 million sheep, goats, 
horses, cattle, camels and yaks in 1990 to about 45 
million animals in 2012 (National Statistical Office of 
Mongolia - NSO, 2012). The increase in animal popula- 
tion was mainly caused by socio-economic changes 
(particularly the breakdown of the Soviet Union), as the 
country's centralized economy turned into a market 
economy (Gong Li et al, 2000). The resulting unemploy- 
ment drove parts of the population back into the country 
side to subsist by increasing the domestic animal herds 
(Reading et al, 2006). 
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Within the scientific community, there has been 
much debate (Vetter, 2005) about whether the observed 
decline in steppe vegetation is related to climate, or 
over-grazing, or both (Fernandez-Gimenez & Allen- 
Diaz, 1999; Opp & Hilbig, 2003). Overgrazing has been 
reported as a major factor of desertification in the 
neighboring province of Inner Mongolia (China) (Gong 
Li et al, 2000), where decline in available grassland 
combined with increase in herd sizes has increased the 
pressure on the remaining pastures, thus accelerating 
desertification (Yiruhan et al, 2001). On the other hand. 
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extreme winter conditions (so called dzud), drought 
and fire, frequently diminish herd sizes, as in 2009/ 
2010 when more than 8 million animals were lost (NSO, 
2012 ). 

The observed decline in grassland vegetation could 
have significant social and ecological impacts region- 
ally and even globally. For instance, related hydrologi- 
cal changes of the land surface have been shown to 
affect regional climate (Gong Li et al., 2000), soil proper- 
ties (Yong-Zhong et al., 2005) and water table (Jia et al., 
2006). The Mongolian steppe ecosystems play an 
important role in mitigating regional and even global 
climate variation through their interaction with the atmo- 
sphere (Yatagai & Yasunari, 1995). In addition, grassland 
ecosystems can act as either a significant sink or source 
for atmospheric C0 2 , thereby potentially contributing to 
or slowing global climate change (Frank, 2002). 

One way to evaluate grassland dynamics at broader 
scales is through time series of remote sensing data. 
Historically, 8 km Normalized Difference Vegetation 
Index (NDVI) data from the National Ocean and Atmo- 
spheric Administration advanced very high resolution 
radiometer (AVHRR) have been used to investigate 
large-scale spatial and temporal patterns of vegetation 
response to climate (Myneni et al., 1997a, b) and photo- 
synthetic capacity (Tucker et al., 2001) since the early 
1980's. AVHRR' s ability to detect subtle changes in 
vegetation dynamics was limited due to instrument 
noise, atmospheric, snow and cloud effects, which 
made it difficult to obtain significant trends in vegeta- 
tion over time (Fensholt & Proud, 2012). Arguably, the 
advent of Moderate Resolution Imaging Spectro- 
radiometer (MODIS) providing near daily coverage at 
250 m-1 km resolution has revolutionized earth system 
science from space, however, some challenges remain 
with respect to removal of atmospheric and cloud 
effects, at least for some regions of the world, including 
tropical (Hilker et al., 2012) and sparsely vegetated 
areas (Steinberg et al., 2006). Recently, a new multi-angle 
implementation of atmospheric correction algorithm 
(MAIAC) (Lyapustin et al., 2011a,b) has shown promise 
to overcome some of these challenges using an advanced 
radiative transfer model for aerosol retrieval, atmo- 
spheric correction and cloud screening (Lyapustin & 
Knyazikhin, 2002; Lyapustin et al., 2012). For instance, in 
previous research across 2 million square kilometers of 
the Amazon basin, Hilker et al. (2012), demonstrated that 
MAIAC can reduce uncertainties in surface reflectance 
up to 10-fold compared to the standard MODIS surface 
reflectance. 

In this study, we use MAIAC to investigate changes 
in vegetation dynamics between 2002 and 2012 based 
on daily MODIS observations. Our objectives are first, 
to derive trends in vegetation greenness between 2002 


and 2012 based on the NDVI (Tucker, 1979) and second, 
to compare these trends to changes in climate variables 
and livestock numbers to investigate potential causes 
for shifts in vegetation density. 

Materials and methods 

Study area 

Our study area encompasses two MODIS tiles (h25v04 and 
h24v04), an area of 2 million square kilometers, spanning 
from roughly 40°N to 50°N in latitude and from 90°E to 
120°E in longitude. Mongolia's climate is semi-arid and 
markedly continental with generally extremely cold, dry 
winters and warm, wet summers. Mean annual temperatures 
range from —1.7 °C in the mountain-steppe to 4.8 °C in the 
desert-steppe, and average annual precipitation ranges 
between 90 mm and 230 mm (Fernandez-Gimenez & Allen- 
Diaz, 1999), most of which falls from fune through Septem- 
ber (Jia et al., 2006). High mountain ranges isolate the coun- 
try from the influence of the Ocean and the Siberian 
anticyclone determines winter climate with its low tempera- 
tures and precipitation rates. Vegetation zones depend on 
altitude, rainfall and soil type and span alpine tundra (3.0% 
of total area), mountain taiga (4.1%), mountain steppe 
(25.1%), steppe (26.1%), desert steppe (27.2%) and desert 
(14.5%) (Hilbig, 1995). Roughly, 124.3 million ha or 79% of 
land area are covered by grassland and about 10% are 
covered by forest or shrub land. 

To characterize changes in climate over the last 30 years, 
we obtained monthly mean values of precipitation and 
temperature from the 64 available weather stations across the 
study area between 1980 and 2012 (National Agency for Mete- 
orology, Hydrology & Environment Monitoring, http://env. 
env.pmis.gov.mn/). The stations are well distributed across 
the country, spanning a range of different climate zones and 
vegetation types (Fig. 1). 

Animal populations 

Official estimates of animal populations were provided on a 
per-province basis by the National Statistical office (NSO) of 
Mongolia. NSO conducts annual surveys of total numbers of 
horses, cattle, sheep, goats and camels using data provided by 
ministries, other state organizations and governors of the 
respective provinces. These livestock estimates are also used 
as official numbers for industry, science and technology in the 
country. 

MAIAC data 

Daily MODIS data from the Aqua platform were collected 
between 2002/07/04 and 2012/12/31 and processed using 
the calibrated and geometrically corrected (Level IB) data. 
We selected MODIS data from the Aqua spacecraft to mini- 
mize uncertainties in vegetation trends from the MODIS 
Terra calibration degradation (Wang et al., 2012). Level IB 
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Fig. 1 Plant functional types derived from the MODIS Land cover product (MCD12Q1). The map presents data from 2002. Over 85% 
of the study area was covered by grassland or sparse vegetation. The black dots symbolize the locations of meteorological stations used 
in this study. 


data were processed using the MAIAC algorithm for cloud 
screening and atmospheric correction (Lyapustin & 
Martonchik, 2011; Lyapustin et al., 2011a,b, 2012). MAIAC 
grids MODIS LIB data (MYD021KM) into 1 km pixels, and 
accumulates measurements of the same surface area from 
different orbits (view geometries) for up to 16 days of 
observations for equatorial and up to 4 days for polar 
regions using a moving window approach (Lyapustin et al., 
2012). The MAIAC cloud mask algorithm composes a 
dynamically updated reference clear-sky image of the 
surface from spatial and time series analyses (Lyapustin et al., 
2008). This allows cloud masking to be based on knowledge 
of reference clear-sky reflectance in addition to spectral and 
thermal thresholds, which can improve cloud detection 
considerably (Hilker et al., 2012). Along with better cloud 
recognition, the MAIAC technique improves aerosol retrieval 
(Lyapustin et al., 2012) and hence the atmospheric correction 
(Lyapustin et al., 2012) by applying advanced radiative trans- 
fer theory (Lyapustin & Knyazikhin, 2001) that does not 
assume a Lambertian reflectance of the surface. The accumu- 
lated multi-angle dataset from multiple overpasses is used to 
derive parameters of the Ross Thick - Li Sparse (RTLS) BRDF 
model (Roujean et al., 1992) for every 1 km grid cell as well as 
the bidirectional reflectance factor (BRF, often called surface 
reflectance) for the last observation. 

Multi-angle implementation of atmospheric correction algo- 
rithm surface products include parameters of BRDF model, 
BRF and albedo at 1 km resolution as well as spectral BRF at 
500 m resolution. The BRF, or surface reflectance, has a signifi- 
cant variability with the view geometry. To mitigate this 
source of variability for trend analysis, we normalized BRFs 
for the fixed view geometry of sun zenith angle 9 S = 45° and 
nadir view (0 V = 0°) using the known BRDF model: 

_ BRF x RTLS(Q s , 6 V , Acj)) 
n RTLS(6 s = 45°, 6 V = 0°, A(/> = 0°) U 

where A(f> is relative azimuth. Such normalization (e.g., 
Lyapustin et al., 2012) was shown to reduce geometric vari- 
ability by a factor of 3-6 thereby helping to isolate seasonal 
changes and multi-year cycles and trends. 


Land cover types 

Different vegetation types have different ecological capacities 
to accommodate stresses induced by biotic or abiotic factors. 
Therefore, it makes sense to investigate changes in vegetation 
cover and their causes separately for different land cover 
types. In this study, we used the MODIS land cover product 
(MCD12Q1, type 5) (Friedl et al., 2010) to stratify time series of 
MAIAC observations into different plant functional types 
(Bonan et al., 2002). The MODIS land cover algorithm is based 
on an ensemble of supervised classification techniques and 
decision tree analysis to yield global land cover types (Friedl 
et al., 2010). Across Mongolia, the product differentiates 
between six different plant functional types (Fig. 1) however; 
the vast majority of the area is dominated by grassland and 
other sparse vegetation. To prevent changes in plant func- 
tional type from affecting vegetation trends (for instance if a 
pixel moves from one land cover class into another during the 
10 year study period), we selected pixels only from areas 
where the plant functional types remained identical between 
2002 and 2012 for analysis. 

Time series fitting 

Time series of satellite data can be used to derive trends in 
vegetation based on the assumption that changes in the 
measured reflectance are driven by phenological changes on 
the ground (Jonsson & Eklundh, 2004). Principally, change 
detection can be done either using simple threshold techniques 
(Sellers et al., 1994) or function fitting (Jonsson & Eklundh, 
2004). While thresholds are easier to implement and resulting 
algorithms are fast, a functional description of changes in 
reflectance has the advantage in allowing seasonal progression 
to be separated from inter-annual progression (Zhu et al., 
2012). Fitting functions through time series data also reduces 
noise levels by smoothing the signal obtained from temporally 
discrete observations (Jonsson & Eklundh, 2004). Different 
algorithms are available and have been applied (Jonsson & 
Eklundh, 2004; Bradley et al., 2007); here we used a simple 
series of harmonic functions (Zhu et al., 2012) based on Fou- 
rier series (Davis & Sampson, 2002). This method developed 
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by Zhu et al. (2012) has been adapted in this study to be used 
for dense time series and across multiple years. One advan- 
tage of Zhu's harmonic model is that it requires only a rela- 
tively small number of clear observations, making predictions 
less sensitive to missing data. At the same time, the model is 
flexible enough to account for intra and inter-annual changes 
in reflectance (Zhu et al., 2012). These properties are particu- 
larly advantageous during the winter months in Mongolia 
where the number of snow free MODIS observations is limited 
and the ground may not be visible for several weeks or even 
months. Our inter-annual change model (i = 2, 3, 4...N) is 
comprised of sines and cosines to represent differences from 
one year to the next (Zhu et al., 2012): 

f(x) = a 0 + (a t cos ^ x'j + hi sin ^ 

flN + lCOS (olH +bN+1 Koir* 

where x, day of year; N, number of years; T, number of days 
per year (365); aO, coefficient to represent the overall surface 
reflectance; 0 Z ,L, coefficients that capture changes of surface 
reflectance from one year to the next; a N+1 ,b N+ i, coefficients 
that capture the bimodal variations of surface reflectance. 

The algorithm allows quantification of changes from oppor- 
tunistically acquired satellite observations and yields continu- 
ous predictions of reflectance or vegetation indices over time. 
Trends in vegetation from one year to the next can be quanti- 
fied in terms of changes in annual mean, maximum and mini- 
mum reflectance and as integral of the modeled reflectance 
over the year. Our harmonic model was fitted to time series of 
NDVI separately for every pixel (1200 x 600 locations); trends 
were derived as inter-annual changes in the fitted curves. 
Specifically, we obtained the mean value for each 365 day 
interval from the fitted curve and derived inter-annual 
changes in NDVI as difference of these annual means. 

Regional differences in factors driving changes in NDVI 
were investigated using a mixed-effects linear regression 
model (Lindstrom & Bates, 1990; Pinheiro & Bates, 1995) to 
assess local effects of changes in climate on NDVI. Mixed 
effects models account for fixed and random effects on the 
variability of a response variable by grouping predictor vari- 
ables into fixed categories (Pinheiro & Bates, 1995). In this 
study, climate anomalies were grouped using long term aver- 
ages (2002-2010) of temperature, precipitation and vegetation 
density (NDVI), respectively, to investigate how climate 
change has affected arid vs. humid vegetation zones and 
sparsely vegetated vs. densely vegetated areas. 



Results 

Figure 2 shows total changes in livestock in Mongolia 
since 1980 (National Statistics Office of Mongolia). The 
figure presents estimated populations of camel, goat, 
sheep, cattle and horses across the country. While num- 
bers in cattle and camels remained relatively constant 
throughout the observation period, a steep increase in 
sheep and, in particular, goat population has been 
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Fig. 2 Changes in livestock populations in Mongolia since 1980 
(source: Statistical Office of Mongolia, NSO). The Figure shows 
a steady increase, particularly in sheep and goat since the early 
1990s. The large increase since 2002 represents recovery from 
the massive die-back caused by the dzud conditions between 
1999-2002. 


recorded since the early 1990s. It is important to note 
the differences in species, as at least in situations where 
animal survival is in competition, goats are likely to 
cause significantly more damage than sheep, which in 
turn will cause more damage than cows (Wilson, 1986). 
The upward trend in number of animals was inter- 
rupted twice. The years 1999-2002 saw three dzud 
conditions in a row, accounting for the loss of about 11 
million animals. In 2010, arctic oscillation events 
resulted in prolonged winter conditions, killing another 
8 million heads of livestock. Despite these calamities, 
total animal populations almost doubled between 2002 
and 2012. Animal populations grew mainly in the wes- 
tern provinces of the country, but increases in livestock 
occurred almost everywhere, with the exception of 
urban areas. Figure 3 shows the spatial distribution of 
cumulative changes in animal populations across the 
study area between 2002 (beginning of the MODIS 
record) and 2012. We show cumulative changes rather 
than simple differences because an increase in number 
of livestock early during the observation period exert 
pressure on the ecosystem longer than if it occurred 
during the last year of our study period. Cumulative 
herd sizes in the western part of the country reached 
up to 5.5 million, while the south eastern region saw 
moderate cumulative increases of up to 0.4 million. 

During the same period, annual mean temperature in 
the country increased by about 1.5 °C (Fig. 4) from 
—0.4 °C in the early 1980s to 1.2 °C in 2010 (Tempera- 
ture data were averaged from annual means using all 
64 available meteorological stations shown in Fig. 1). 
Meteorological data also showed a decrease in precipi- 
tation of about 5% from, on average, 204 mm yr -1 in 
the 1980s to about 195 mm yr -1 in the early and mid 
2000s. Most of the warming occurred in the northern 
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Fig. 3 Cumulative changes in herd sizes (in million animals) since 2002. 



Fig. 4 Trends in average temperature and total annual precipi- 
tation across Mongolia between 1980-2010 (mean of all available 
weather stations). 


and central provinces of the country (Fig. 5), with mean 
annual increases in temperature of up to 0.1 °C. Central 
Mongolia also experienced the largest decrease in 
annual precipitation, on average 0.4 mm yr _1 over the 
30 year period. 

Figure 6 illustrates an example of a harmonic model 
used to fit NDVI series. This example originates from 
the northern part of the study area. Seasonal changes 
were large for this particular example, reflecting 


Mongolia's continental climate. NDVI ranged from about 
0.7 at the height of the growing season to 0.2 and the end 
of the growing season. The Figure also shows gaps in the 
analysis in the winter months, due to snow cover. 

The quality of the MAIAC time series of observation 
is discussed in detail in Appendix SI. The harmonic 
model chosen in this study fitted the time series in this 
example well, although it somewhat underestimated 
the peak of the growing season. Figure 7 shows the 
standard error (SE) of the harmonic models fitted at 
each pixel location across the study area. The harmonic 
functions fitted the MODIS observations well. The SE 
was highest in the northern part for the study area 
(0.04-0.06) and lowest in the southern part of the coun- 
try with sparse vegetation cover (0.01). Figure 8 shows 
changes in yearly mean NDVI (as linear trend between 
2002 and 2012), derived using the annual mean of the 
harmonic model. Figure 8(a) presents annual mean 
NDVI as observed 2002. Mean NDVI was highest in the 
northern part of the study area reaching up to 0.5, 
whereas the southern region bordering the Gobi desert 
showed much lower values (0.1-0.2). Figure 8(b) 
presents changes in NDVI (delta NDVI) observed 
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Fig. 5 Trends in temperature and precipitation between 1980 and 2010 by province. Changes were estimated as slope of linear trends 
for each weather station. 
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Fig. 6 Example of a harmonic fit using combined MODIS Aqua NDVI data acquired between 2000 and 2012. The gray dots represent 
NDVI MODIS MAIAC observations of NDVI, the black line is the fitted model between 2002 and 2012. 


I WE 96 WE 100 WE 105 WE HOWE 115 WE 





45 WN- 



•50WN 


-45WN 


-40 WN 


HJtT£ 95 WE 1D0 a 0'0’ , E 105 WE HOWE 11 5 WE 


Fig. 7 Standard error (SE) of the harmonic model for BRF normalized NDVI data. The SE was around 5% in the northern part and 
1-2% in the southern part of the study area. 
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Fig. 8 (a) Mean NDVI over the 2002 growing season (May-September), (b) Changes in mean NDVI between 2002 and 2012, derived 
using the harmonic model. 


between 2002 and 2012. While NDVI was relatively 
stable or increased in some areas in the north of the 
country, a decrease of up to 0.05 in yearly mean NDVI 


occurred over much of the southern provinces, espe- 
cially in the most southern part of the country. The 
decline was especially pronounced across the southern 
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edge of the grassland (compare Fig. 1) and along the 
transition into barren /sparse vegetation functional 
type. This region experienced 25-50% decrease in 
NDVI compared to the 2002 values (Fig. 8a). Across the 
whole study area, NDVI decreased by 0.025 which cor- 
responds to a decrease of about 12% of total mean 
NDVI observed in 2002. 

Changes in NDVI were strongly related to cumula- 
tive increases in animal population. Figure 9 shows the 
relationship between mean annual changes in NDVI 
between 2002 and 2012 and cumulative change in herd 
size as reported for every province. Figure 9(a) shows 
the decline in grassland vegetation (plant functional 
type). Figure 9(b) shows the respective changes for the 
barren/sparse vegetation type. A strong linear relation- 
ship between cumulative increase in animal population 


and decrease in NDVI was found for both plant func- 
tional types, with r 2 values of 0.85 and 0.74 for grass- 
land and sparse vegetation /barren types, respectively. 
Note that some of the very northern and western prov- 
inces were not fully contained in the two MODIS tiles 
used in this analysis. NDVI decreased by up to 0.07 in 
some provinces with cumulative increase in herd sizes 
were up to 5 million heads of livestock. For the sparse 
vegetation type, the corresponding decline in NDVI 
was about 0.02, presumably because those areas 
already had less vegetation cover at the start of the 
observation period. 

Figure 10 shows the relationship between changes in 
NDVI and precipitation (a) and changes in NDVI and 
temperature, (b) Data are reported as deviation from 
the average annual mean between 2002 and 2012. NDVI 




Fig. 9 Relationship between cumulative increase in animal population and decline in NDVI. 


Precipitation Temperature 




Fig. 10 (a) Relationship between changes in NDVI (as deviation from the long term average between 2002 and 2012) and changes in 
total summer precipitation (here defined as precipitation between May and September). Figure 10 (b) Relationship between changes in 
NDVI (as deviation from the long term average between 2002 and 2012) and deviation in mean annual temperature from the long term 
mean (as deviation from the long term average between 2002 and 2012). 
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values varied by about 0.03 around the long term mean 
(2002-2012), whereas precipitation and temperature 
deviated by about ±70 mm and ±2 °C, from their long 
term (2002-2012) mean values, respectively. As a 
whole, a moderate positive relationship was found 
between increase in summer precipitation (May 
through September) and NDVI (r 2 = 0.30; P < 0.05) and 
a weak negative relation was found between increase 
in temperature and NDVI (r 2 = 0.11; P < 0.05). The 
values shown here were averaged by provinces to make 
the results more comparable to those presented in 
Fig. 9. 

While climatic changes explained about 30% of vege- 
tation decline across the country as a whole, our mixed 
effects model showed considerable regional differences 
in the climatic effects on NDVI (Fig. 11). From the fixed 
effects that were investigated in this study (average 
vegetation density, average precipitation rate and aver- 
age temperature), only vegetation density had a signifi- 
cant impact on the relationship between climate 
anomalies and NDVI anomalies. In densely vegetated 
areas (mean NDVI > 0.60), anomalies in precipitation 
explained about to 50% of the anomalies in NDVI 
(P < 0.05). In sparsely vegetated areas (NDVI < 0.2) 
less than 10% of changes in NDVI were explained by 
changes in precipitation (Fig. 11(a), P < 0.05). Similarly, 
temperature anomalies explained about 20% of changes 
in NDVI (P < 0.05) in the densely vegetated north and 
north-east of the country, while the effect on most of 
the grassland area in the south and west was insignifi- 
cant. We also related NDVI anomalies to changes in 
precipitation (Fig. 11a) and temperature (Fig. lib), 
separately for each meteorological station. For most 
meteorological stations located in the south, only weak 
relationships were found between changes in NDVI 
and changes in precipitation and temperature, while 
climate anomalies generally explained more of the 


inter-annual variability in NDVI at meteorological 
stations located in denser vegetated areas (Fig. 11). 

Discussion 

We have observed and analyzed changes in grassland 
vegetation across Mongolia using MODIS satellite data 
between 2002 and 2012. A key finding of this study is 
that the degradiation in grasslands over the last 10 
years across Mongolia (on average by about 12% of 
mean NDVI observed in 2002) is mainly related to an 
increase in domestic animal populations, and to a more 
limited extent, to changes in precipation patterns. 
While it is difficult from remote sensing data alone to 
assess to what extend this decrease in NDVI translates 
into decrease in plant biomass, it can be assumed that 
the decline should be in about the same range, at least 
in the southern regions, as the relationship between 
NDVI and vegetation leaf area is linear for leaf areas of 
less than about 4 m 2 m~ 2 (Carlson & Ripley, 1997; 
Myneni et al., 1997b; Turner et al., 1999). The observed 
degradiation in vegetation is well supported by previ- 
ous studies, based mainly on data acquired from 
ground observations (Fernandez-Gimenez & Allen-Diaz, 
1999; Stumpp et al., 2005), but also AVHRR satellite 
observations (Sternberg et al., 2011; Liu et al., 2013). We 
found most degradiation in the central and southern 
part of the country, particularly in the transition zone 
between grassland and barren /sparse vegetation with 
up 40% decline compared to 2002. This finding suggests 
that grasslands are retreating rapidly in the transition 
zone, thus leading to expansion of deserted areas 
(Sternberg et al., 2011). Changes in summer precipita- 
tion seem to play some role in this decline (Fig. 10), 
especially in the central and eastern part of the country 
(Figs 5 and 11). While significant regional differences 
were found on the impact of changes in precipitation 


■SOWN 


) H N 


■40 D 00"N 

Fig. 11 Regional differences in the variance explained in NDVI across the country using mixed effects models. The map colors indicate 
the coefficient of determination (r 2 ) for the relationship between precipitation anomalies and NDVI anomalies (Fig. 11a) and tempera- 
ture anomalies and NDVI anomalies (Fig. lib). In both cases long term mean NDVI was used as a grouping variable (steps of 0.2 from 
0 to 1). The color of the dots represents the strength of the relationship (r 2 ) between precipitation anomalies and NDVI anomalies 
(Fig. 11a) and temperature anomalies and NDVI anomalies (Fig. lib) at individual meteorological stations. 
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patterns on vegetation decline, our findings suggest 
that over-grazing is the single most important reason 
for desertification of the Mongolian Steppe. This con- 
clusion is supported by similar studies from neighbor- 
ing Inner Mongolia (China), that have identified over- 
grazing as a major cause of desertification (Gong Li 
et al, 2000). Zhao et al. (2005) showed that heavy graz- 
ing results in a considerable decrease in vegetation 
cover, root biomass, and an increase in animal hoof 
impacts. Batjargal (1997) identified livestock grazing as 
most prevalent human-related activity contributing to 
land degradation in the area. Liu et al. (2013) reported a 
steady decline in vegetation across Mongolia between 
1998 and 2008. They attributed changes mostly to cli- 
matic shifts, based on radar derived vegetation optical 
depth. While their results seemingly contradict the 
findings of this work, the discrepancies may be 
explained by regional differences in the factors driving 
decline in vegetation and the complex interactions 
between them. For instance, Liu et al. (2013) considered 
only increases in goat populations, which were mostly 
limited to central Mongolia (Liu et al., 2013). In this 
region, our study has shown the strongest decline in 
precipitation (Fig. 5), which had a clear negative effect 
on vegetation growth (Fig. 11). We found similar 
results for the north-eastern part of Mongolia, changes 
in precipitation, however, could not explain the vegeta- 
tion decline observed across large areas in the southern 
and western part of the country (Figs 8 and 11). Over- 
grazing and climate variabity therefore need to be 
understood as two extremes along a gradient, rather 
than contradicting theories (Illius & O'connor, 1999). 
For instance, our results suggest that degradation in 
central Mongolia is largely caused by overgrazing, 
however, changes in summer preciptation have 
occurred (Fig. 5) and they had a clear impact on regen- 
eration of grasslands (Fig. 11). Such changes may 
become more important once critical thresholds are 
reached and exceeded. 

Our analysis benefited from using an improved sur- 
face reflectance product (MAIAC), thereby confirming 
the value of previous work (Hilker et al., 2012) by 
demonstrating that high quality reflectance data can 
be consistently derived from this product (Figs 6 and 7). 
While 11 years of time series from MODIS/ MAIAC is 
relatively short, results presented in the appendix 
have shown that MAIAC processed MODIS data 
allow a more accurate assessment of land degradation 
compared to for instance AVHRR. In addition, results 
presented in Fig. 2 show that increase in livestock has 
not been linear over the last 30 years; in fact, the vast 
majority of it has occurred since the early 2000s. As a 
result, it can be assumed that the majority of live- 
stock-related disturbance has happened since then. 


The higher error in MAIAC derived NDVI in the 
northern part of the country could be due to snow 
cover, which led to extended periods during which no 
valid estimates of ground cover could be made), or 
the inability of multi-harmonic model to accurately fit 
the seasonality (e.g., underestimation of NDVI peaks 
in Fig. 6). 

The selected time series model fitted the observations 
generally well, as shown in the example of Fig. 6 and 
by the standard error presented in Fig. 7. These results 
compare well also to previous studies based on Landsat 
observations (Zhu et al., 2012). While the approach 
applied in this study was successful, it may be limited 
in other areas, as it makes the assumption that vegeta- 
tion has only a one season cycle, and it excluded places 
where a change in vegetation type may have occurred. 
While several other techniques are available (Jonsson & 
Eklundh, 2004; Huang & Wu, 2008) that do not apply 
such assumptions and exclusions, our technique is 
more likely to be insensitive to data gaps, for instance 
as a result of prolonged snow cover, than alternative 
procedures. 

While we were able to associate recent decline in 
grassland with grazing, it is unknown to what extend 
changes in climate may have contributed to predispos- 
ing the steppe ecosystem to stress. Temperatures have 
increased considerably over the last 30 years across the 
country, but its effect on vegetation was limited, at least 
over the last 10 years. It is reasonable to assume that in 
semi-arid environments, water, rather than tempera- 
ture limit growth (see Fig. 10a). An increase in temper- 
ature could therefore affect plant growth in a negative 
way, because it will enhance in latent heat flux, there- 
fore increasing the amount of water evaporating from 
the surface. 

We have shown a clear connection between increases 
in animal population and decline in vegetation density. 
This trend has continued despite two major die-offs of 
animal populations during harsh winter conditions. We 
recognize that our study was limited to studying very 
recent climate related impacts on grassland vegetation. 
Longer-term trends in climate were not evaluated and 
may play a significant role, particularly if extreme 
events are becoming more frequent. 

Acknowledgements 

We would like to thank Dr. Paul Doescher (Oregon State 
University) for providing funding for Ms. Natsagdorj to work on 
this project. RHW contributions to this paper are an extension of 
his research supported by NASA grant NNX11A029G. The work 
of Drs. Lyapustin and Wang was supported by the NASA 
Science of Terra and Aqua Program. Thanks to the National Sta- 
tistical Office of Mongolia and the Institute of Metrological 
Hydrology for providing meteorological and animal population 


© 2013 John Wiley & Sons Ltd, Global Change Biology , 20 , 418-428 


OVERGRAZING CAUSES MONGOLIA GRASSLAND DECLINE 427 


data used in this study. We would like to thank Drs. Woodcock 
and Zhu for helpful discusions regarding the harmonic func- 
tions used to fit the time series of satellite data in this study. 

References 

Batjargal Z (1997). In Proceedings of the International Workshop on Rangeland 
Desertification. RALA Report (Vol. 200). 

Bonan GB, Levis S, Kergoat L, Oleson KW (2002) Landscapes as patches of plant func- 
tional types: an integrating concept for climate and ecosystem models. Global 
Biogeochemical Cycles, 16, 5-1-5-23. 

Bradley BA, Jacob RW, Hermance JF, Mustard JF (2007) A curve fitting procedure to 
derive inter-annual phenologies from time series of noisy satellite NDVI data. 
Remote Sensing of Environment, 106, 137-145. 

Carlson TN, Ripley DA (1997) On the relation between NDVI, fractional vegetation 
cover, and leaf area index. Remote Sensing of Environment, 62, 241-252. 

Davis J, Sampson R (2002) Statistics and Data Analysis in Geology. (3rd edn). Wiley, 
New York. 

Fensholt R, Proud SR (2012) Evaluation of Earth observation based global long term 
vegetation trends-Comparing GIMMS and MODIS global NDVI time series. 
Remote sensing of Environment, 119, 131-147. 

Fernandez-Gimenez ME, Allen-Diaz B (1999) Testing a non-equilibrium model of 
rangeland vegetation dynamics in Mongolia. Journal of Applied Ecology, 36, 871- 
885. 

Fontana FMA, Coops NC, Khlopenkov KV, Trishchenko AP, Riffler M, Wulder MA 
(2012) Generation of a novel 1 km NDVI data set over Canada, the northern 
United States, and Greenland based on historical AVHRR data. Remote Sensing of 
Environment, 121, 171-185. 

Frank A (2002) Carbon dioxide fluxes over a grazed prairie and seeded pasture in the 
Northern Great Plains. Environmental Pollution, 116, 397-403. 

Frey RA, Ackerman SA, Liu Y, Strabala KI, Zhang H, Key JR, Wang X (2008) Cloud 
detection with MODIS. Part I: improvements in the MODIS cloud mask for collec- 
tion 5. Journal of Atmospheric and Oceanic Technology, 25, 1057-1072. 

Friedl MA, Sulla-Menashe D, Tan B, Schneider A, Ramankutty N, Sibley A, 
Huang X (2010) MODIS Collection 5 global land cover: algorithm refinements 
and characterization of new datasets. Remote Sensing of Environment, 114, 168- 
182. 

Gong Li S, Harazono Y, Oikawa T, Zhao HL, Ying He Z, Chang XL (2000) Grassland 
desertification by grazing and the resulting micrometeorological changes in Inner 
Mongolia. Agricultural and Forest Meteorology, 102, 125-137. 

Hilbig W (1995) The Vegetation of Mongolia. SPB Academic Publishing, Amsterdam. 

Hilker T, Lyapustin AI, Tucker CJ, Sellers PJ, Hall FG, Wang Y (2012) Remote sensing 
of tropical ecosystems: atmospheric correction and cloud masking matter. Remote 
Sensing of Environment, 127, 370-384. 

Huang NE, Wu Z (2008) A review on Hilbert-Huang transform: method and its appli- 
cations to geophysical studies. Reviews of Geophysics, 46, RG2006. 

Illius AW, O'connor TG (1999) On the relevance of nonequilibrium concepts to arid 
and semiarid grazing systems. Ecological Applications, 9, 798-813. 

Jia B, Zhou G, Wang Y, Wang F, Wang X (2006) Effects of temperature and soil water- 
content on soil respiration of grazed and ungrazed Leymus chinensis steppes. 
Inner Mongolia. Journal of Arid Environments, 67, 60-76. 

Jonsson P, Eklundh L (2004) TIMESAT-a program for analyzing time-series of satel- 
lite sensor data. Computers & Geosciences, 30, 833-845. 

Kalinina A (1974) Main Types of Pastures of Mongolia, their Structure and Productivity. 
Biological Resources and Natural Conditions of the Mongolian People's Republic. IE 
Izdatalestvo Nauka, Leningrad, [in Russian]. 

Kawamura K, Akiyama T, Yokota H, et al. (2007) Monitoring of forage conditions 
with MODIS imagery in the Xilingol steppe. Inner Mongolia. International Journal 
of Remote Sensing, 26, 1423-1436. 

Lavrenko EM, Karamysheva ZV (1993) Steppes of the Former Soviet Union and Mongolia. 
Elsevier, Amsterdam. 

Li X, Strahler AH (1986) Geometric-optical bidirectional reflectance modeling of a 
conifer forest canopy. IEEE Transactions on Geoscience and Remote Sensing, 6, 
906-919. 

Li SG, Asanuma J, Eugster W, Kotani A, Liu JJ, Urano T, Sugita M (2005) Net ecosys- 
tem carbon dioxide exchange over grazed steppe in central Mongolia. Global 
Change Biology, 11, 1941-1955. 

Lindstrom M, Bates D (1990) Nonlinear mixed effects models for repeated measures 
data. Biometrics, 46, 673-687. 


Liu YY, Evans JP, McCabe MF, de Jeu RAM, van Dijk AIJM, Dolman AJ, Saizen I 
(2013) Changing climate and overgrazing are decimating Mongolian steppes. PloS 
one, 8, e57599. 

Lyapustin A, Knyazikhin Y (2001) Green's function method for the radiative transfer 
problem. I. Homogeneous non-Lambertian surface. Applied optics, 40, 3495-3501. 

Lyapustin A, Knyazikhin Y (2002) Green's function method in the radiative transfer 
problem. II. Spatially heterogeneous anisotropic surface. Applied optics, 41, 5600- 
5606. 

Lyapustin A, Wang Y, Frey R (2008). An automatic cloud mask algorithm based on 
time series of MODIS measurements. Journal of Geophysical Research, 113, D16207. 

Lyapustin A, Martonchik J, Wang Y, Laszlo I, Korkin S (2011a) Multiangle implemen- 
tation of atmospheric correction (MAIAC): 1. Radiative transfer basis and look-up 
tables. Journal of Geophysical Research, 116, D03210. 

Lyapustin A, Wang Y, Laszlo I, Kahn R, Korkin S, Remer L, Reid JS (2011b). Multian- 
gle implementation of atmospheric correction (MAIAC): 2. Aerosol algorithm. 
Journal of Geophysical Research, 116, D03211. 

Lyapustin A, Wang Y, Laszlo I, Hilker T (2012) Multi-angle implementation of atmo- 
spheric correction for MODIS (MAIAC). Part 3: atmospheric correction. Remote 
Sensing of Environment, 127, 385-393. 

Myneni RB, Keeling CD, Tucker CJ, Asrar G, Nemani RR (1997a) Increased plant 
growth in the northern high latitudes from 1981 to 1991. Nature, 386, 698-702. 

Myneni RB, Ramakrishna R, Nemani R, Running SW (1997b) Estimation of global leaf 
area index and absorbed PAR using radiative transfer models. IEEE Transactions 
on Geoscience and Remote Sensing, 35 (6), 1380-1393. 

National Statistical Office of Mongolia (2012) Available at: http://web.nso.mn:8080/ 
userdata/Dialog/varval.asp?ma=MAAl&ti=&path=../Database/Mongolian/ 
MAA/&lang=l (accessed 4 April 2013). 

Opp C, Hilbig W (2003) The impact of overgrazing on natural pastures in Mongolia 
and Tyva. Berliner Paldobiologische Abhandlungen, 2, 96e98. 

Pacyna A (1986) Vegetation of the Sant valley in the Khangai mountains (Mongolia). 
Fragmenta Floristica Geobotanica, 30, 313-451. 

Pinheiro JC, Bates DM (1995) Approximations to the log-likelihood function in the 
nonlinear mixed-effects model. Journal of Computational and Graphical Statistics, 4, 
12-35. 

Pinzon J, Brown ME, Tucker CJ (2005) Satellite time series correction of orbital drift 
artifacts using empirical mode decomposition. In: Hilbert-Huang Transform: Intro- 
duction and Applications, (ed. Huang N), pp. 167-186. World Scientific, New Jersey. 

Reading RP, Bedunah DJ, Amgalanbaatar S (2006) Rangelands of Central Asia: Proceed- 
ings of the Conference on Transformations, Issues, and Future Challenges, In: USDA 
Forest Service Proceedings. Fort Collins, CO. 

Roujean J-L, Leroy M, Deschamps P-Y (1992) A bidirectional reflectance model of the 
Earth's surface for the correction of remote sensing data. Journal of Geophysical 
Research, 97, 20455. 

Ross I (1981) The Radiation Regime and Architecture of Plant Stands. Dr. W. Junk, Nor- 
well, MA 

Sellers PJ, Tucker CJ, Collatz GJ, Los SO, Justice CO, Dazlich DA, Randall DA (1994) 
A global 1° by 1° NDVI data set for climate studies. Part 2: The generation of glo- 
bal fields of terrestrial biophysical parameters from the NDVI. International Journal 
of Remote Sensing, 15, 3519-3545. 

Steinberg DC, Goetz SJ, Hyer EJ (2006) Validation of MODIS F/sub PAR/ products in 
boreal forests of Alaska. IEEE Transactions on Geoscience and Remote Sensing, 44, 
1818-1828. 

Sternberg T, Tsolmon R, Middleton N, Thomas D (2011) Tracking desertification on 
the Mongolian steppe through NDVI and field-survey data. International Journal of 
Digital Earth, 4, 50-64. 

Stumpp M, Wesche K, Retzer V, Miehe G (2005) Impact of grazing livestock and dis- 
tance from water source on soil fertility in southern Mongolia. Mountain Research 
and Development, 25, 244-251. 

Sugita M, Asanuma J, Tsujimura M, Mariko S, Lu M, Kimura F, .. & Adyasuren, T. 
(2007) An overview of the rangelands atmosphere-hydrosphere-biosphere inter- 
action study experiment in northeastern Asia (RAISE). Journal of Hydrology, 333, 
3-20. 

Tucker CJ (1979) Red and photographic infrared linear combinations for monitoring 
vegetation. Remote Sensing of Environment, 8, 127-150. 

Tucker CJ, Slayback DA, Pinzon JE, Los SO, Myneni RB, Taylor MG (2001) Higher 
northern latitude normalized difference vegetation index and growing season 
trends from 1982 to 1999. International Journal of Biometeorology, 45, 184-190. 

Tucker CJ, Pinzon JE, Brown ME (2004) Global Inventory Modeling and Mapping Studies, 
NA94aprl5b.nll-VIg. Global Land Cover Facility, University of Maryland, College 
Park, MD. 


2013 John Wiley & Sons Ltd, Global Change Biology , 20 , 418-428 


428 T. HILKER et ah 


Tucker C, Pinzon J, Brown M et al. (2005) An extended AVHRR 8-km NDVI dataset 
compatible with MODIS and SPOT vegetation NDVI data. International Journal of 
Remote Sensing, 26 , 4485-4498. 

Turner DP, Cohen WB, Kennedy RE, Fassnacht KS, Briggs JM (1999) Relationships 
between leaf area index and Landsat TM spectral vegetation indices across three 
temperate zone sites. Remote Sensing of Environment, 70, 52-68. 

UNEP (2002) Mongolia: State of the environment. United Nations Environment Pro- 
gramme Pathumthani, Thailand. 

Vermote EF, Kotchenova SY, Ray JP (2008) MODIS Surface Reflectance User's Guide 
v 1.2. Available at: http://modis-sr.ltdri.org (accessed 27 March 2013). 

Vetter S (2005) Rangelands at equilibrium and non-equilibrium: recent developments 
in the debate. Journal of Arid Environments, 62 , 321-341. 

Wallis de Vries MF, Manibazar N, Diigerlham S (1996) The vegetation of the forest- 
steppe region of Hustain Nuruu, Mongolia. Vegetatio, 122, 111-127. 

Wang D, Morton D, Masek J, Wu A, Nagol J, Xiong X, Wolfe R (2012) Impact of sen- 
sor degradation on the MODIS NDVI time series. Remote Sensing of Environment, 
119 , 55-61. 

Wilson AD (1986) Principles of Grazing Management Systems, Proceedings of the Inter- 
national Rangeland Congress. Adelaide, Australia, pp. 221-225. 

Yatagai A, Yasunari T (1995) Interannual variations of summer precipitation in the 
arid /semi-arid regions in China and Mongolia: their regionality and relation to 
the asian summer monsoon. Journal of the Meteorological Society of Japan, 73, 
909-923. 

Yiruhan I, Hayashi I, Nakamura T, Shiyomi M (2001) Changes in floristic composition 
of grasslands according to grazing intensity in inner Mongolia, China. Journal of 
Japanese Society of Grassland Science, 47, 362-369. 


Ykhanbai H, Bulgan E, Beket U, Vernooy R, Graham J (2004) Reversing grassland 
degradation and improving Herders' livelihoods in the Altai Mountains of 
Mongolia. Mountain Research and Development, 24 , 96-100. 

Yong-Zhong S, Yu-Lin L, Jian-Yuan C, Wen-Zhi Z (2005) Influences of continuous 
grazing and livestock exclusion on soil properties in a degraded sandy grassland. 
Inner Mongolia, northern China. Catena, 59 , 267-278. 

Zhao H-L, Zhao X-Y, Zhou R-L, Zhang T-H, Drake S (2005) Desertification processes 
due to heavy grazing in sandy rangeland. Inner Mongolia. Journal of Arid Environ- 
ments, 62 , 309-319. 

Zhu Z, Woodcock CE, Olofsson P (2012) Continuous monitoring of forest disturbance 
using all available Landsat imagery. Remote Sensing of Environment, 122 , 75-91. 


Supporting Information 

Additional Supporting Information may be found in the 
online version of this article: 

Appendix SI. Comparison between MAIAC and conven- 
tional NDVI. 


© 2013 John Wiley & Sons Ltd, Global Change Biology , 20 , 418-428 


